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We present version 2 of QuTiP, the Quantum Toolbox in Python. Compared to the preceding version [Comput. Phys. 
Comm. 183 (2012) 1760], we have introduced numerous new features, enhanced performance, made changes in the 
Application Programming Interface (API) for improved functionality and consistency within the package, as well as 
increased compatibility with existing conventions used in other scientific software packages for Python. The most sig- 
nificant new features include efficient solvers for arbitrary time-dependent Hamiltonians and collapse operators, support 
for the Floquet formalism, and new solvers for Bloch-Rcdfield and Floquet-Markov master equations. Here we introduce 
these new features, demonstrate their use, and give a summary of the important backward-incompatible API changes 
introduced in this version. 
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Restrictions: Problems must meet the criteria for using the 
master equation in Lindblad, Floquet-Markov, or Bloch- 
Redfield form. 

Running time: A few seconds up to several tens of hours, 
depending on size of the underlying Hilbert space. 



1. Introduction 

The Quantum Toolbox in Python (QuTiP), is a generic 
framework for numerical simulation and computation of 
the dynamics of both open and closed quantum systems. 
This framework follows an object-oriented design that 
makes the programming of a quantum mechanical problem 
an intuitive process that closely follows the corresponding 
mathematical formulation. Although implemented in the 
interpreted programming language Python, the use of the 
NumPy and SciPy scientific libraries, and selective opti- 
mization using Cython, allows QuTiP to achieve perfor- 
mance that matches, or in many cases, exceeds that of 
natively compiled alternatives. At the same time, QuTiP 
provides a significantly more convenient and flexible pro- 
gramming environment that is easy to learn, and is well 
suited for use in the classroom. 

Since the first major-version release of QuTiP [1] (the 
1.x series), active development of the framework has re- 
sulted in a significant number of new features and perfor- 
mance enhancements, that have culminated in the second 
major-version release. Here we briefly describe the most 
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significant changes and additional functionality introduced 
in this most recent release. 

This paper is organized as follows. In Sec. [2] we high- 
light the important API changes introduced in going from 
QuTiP 1.x to versions 2.x and higher. Section[3]details the 
primary new features included in this latest version. To 
illustrate the new functionality in QuTiP, Sec. [4] contains 
a selection of examples that highlight how these functions 
are used in numerical quantum simulations. Finally, a list 
of all new user-accessible functions, including a brief de- 
scription, is given in Appendix A. 

2. API changes 

Here we list the backward-incompatible changes in the 
API of QuTiP 2 as compared to the previous version 
(1.1.4) described in Ref. pQ. These changes are impor- 
tant when porting applications and simulations that are 
developed for QuTiP 1.1.4 to QuTiP 2.0 and higher. For 
newly developed simulations, we recommend following the 
documentation and examples for QuTiP 2.1 [2J, in which 
case the following API changes are not relevant. Note that 
this article covers version 2.1.0 of the QuTiP framework, 
and incorporates several features added since the initial 
2.0 release. 

1. All quantum dynamics solvers (mcsolve, mesolve, 
odesolve, essolve, brmesolve, and fmmesolve) now 
return an Odedata instance, that contains all informa- 
tion about the solution (as opposed to data lists or 
Qobjs lists as in QuTiP version 1.1.4). A typical call 
to a time-evolution solver using the new API is: 

sol = solver(H, psiO, tlist, c_ops, e_ops) 

where the sol. expect or sol. states contain the lists 
of expectation values or Qobj instances, respectively, 
that would be returned by the same solver in QuTiP 
1.1.4. With this new API, each solver can optionally 
store additional information in the return object such 
as, for example, the collapse times calculated in the 
Monte Carlo solver. 

2. The name of the function for Lindblad master equation 
solver has been changed from odesolve to mesolve. 
The odesolve function can still be called, however it 
is officially deprecated, and will be removed in a future 
release. Being a QuTiP version 1.x function, odesolve 
does not return an Odedata object. 

3. The order of the return values of the method 
Qobj . eigenstates have been swapped, so that the 
eigenenergies and eigenstates of a Qobj instance op are 
now returned in the following order: 

eigvals, eigkets = op . eigenstates () 

4. Functions for calculating correlations using different 
solvers have now been consolidated under the functions 
correlation and correlation_ss, for transient and 
steady state correlations respectively. Here, the selec- 
tion of the underlying dynamics solver now is specified 



using the optional keyword argument solve that de- 
faults to the Lindblad master equation (mesolve) if it 
is not explicitly specified. For example: 

corr_mat = correlation(H, rhoO, tlist, taulist, 
c_op, A, B, solver="me") 

where solver can be "me", "es" or "mc". 

3. New features 

QuTiP 2 includes a wide variety of new computational 
functions, as well as utility functions for better handling 
of data. Here we give a brief description the new major 
features in QuTiP 2.1. For full documentation of these 
new features, as well as the rest of the QuTiP package, see 
the QuTiP 2.1 Documentation [2|. Examples illustrating 
the usage of these functions can be found in Sec. [3] 

• Support for time-dependent collapse opera- 
tors: We have created a new system for represent- 
ing time-dependent quantum operators used in defin- 
ing system Hamiltonians and collapse operators for 
the Lindblad master equation and the Monte Carlo 
solvers. This allows support for arbitrary time- 
dependent collapse operators (the new method is still 
backwards compatible with the python function call- 
back method used for time-dependent Hamiltonians 
in QuTiP version 1.1.4). This new method of defin- 
ing arbitrary time-dependencies is both more efficient 
and more flexible, allowing for high performance simu- 
lations of arbitrary time-dependent quantum systems. 
In particular, many problems of interest may be com- 
piled at runtime into C code via Cython [jj]. This 
particular feature, and its implementation, will be dis- 
cussed elsewhere [4]. 

• Floquet formalism, Floquet-Markov master 
equation: For periodic time-dependent systems, the 
Floquet formalism can be a useful technique where the 
original time-dependent problem is transformed into a 
time-independent problem using the time-dependent 
Floquet-modes as the basis set. In QuTiP 2.0 we 
added a new module for the Floquet-related decompo- 
sition of time-dependent problems, and the evolution 
of unitary and dissipative dynamics using equations 
of motion and master equations in the Floquet for- 
malism. 

• Bloch-Redfield master equation solver: A new 

quantum dynamics solver for the time evolution ac- 
cording to the Bloch-Redfield master equation is now 
included in QuTiP. While not as efficient as the Lind- 
blad master equation solver, in situations where the 
environment is expressed in terms of its noise power 
spectrum, rather than phenomenological decay and 
dephasing rates used in the Lindblad formalism, the 
Bloch-Redfield master equation has significant advan- 
tages. 
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• Quantum Process Tomography: Quantum pro- 
cess tomography (QPT) [5] is a useful technique for 
characterizing experimental implementations of quan- 
tum gates involving a small number of qubits. It can 
also be a useful theoretical tool that gives insight into 
how a given process transforms density matrices, and 
can be used, for example, to study how noise or other 
imperfections deteriorate quantum gate operations. 
Unlike the fidelity or trace distance, that give a single 
number indicating how far from ideal a gate is, quan- 
tum process tomography gives detailed information as 
to exactly what kind of errors various imperfections 
and losses introduce. 

• Functions for generating random states and 
matrices: It is now possible to generate random kets, 
density matrices, Hamiltonians, and Unitary opera- 
tors. This includes the ability to set the sparsity (den- 
sity) of the resultant quantum object. 

• Support for sparse eigensolvers: The quan- 
tum object (Qobj) methods eigenstates, 
eigenenergies, and groundstate can now use 

sparse eigensolvers for systems with large Hilbert 
spaces. This option is not enabled by default, and 
must be set with the keyword argument sparse. 

• New entropy and entanglement functions: 

Functions for calculating the concurrence, mutual in- 
formation, and conditional entropy have been added. 

• New operator norms: When calculating opera- 
tor norms, one can now select between the following 
norms: trace, Frobius, one, and max. The trace norm 
is chosen by default. For ket and bra vectors, only the 
L2-norm is available. 

• Saving and loading data: Saving and loading of 
quantum objects and array data is now internally 
supported by QuTiP. The storage and retrieval of all 
quantum objects and Odedata objects can be accom- 
plished via the qsave and qload functions, respec- 
tively. In order to facilitate the export of QuTiP 
data to other programs, the f ile_data_store and 
f ile_data_read allow the user to read and write ar- 
ray data, both real and complex, into text files with 
a wide variety of formatting options. 

• Performance improvements: In QuTiP 2.1, nu- 
merous performance optimizations have been imple- 
mented, including more efficient quantum object cre- 
ation, significantly faster ptrace implementation, and 
an improved steadystate solver. 

• Unit tests for verification of installation: The 

installation of QuTiP 2.1 comes with a set of unit 
tests that can be used to verify that the installation 
was successful, and that the underlying routines are 
functioning as expected. 



4. Example scripts featuring new functionality 

In this section we highlight, via examples, several of the 
main features added in QuTiP 2.1 and listed in Sec. [3| 
Although we will demonstrate the use of the new time- 
dependent evolution framework, a full discussion of this 
feature is presented elsewhere [3] . The examples listed be- 
low, as well as a growing collection of additional demon- 
strations, can be found on the QuTiP website [2] , or run af- 
ter installing QuTiP using the demos function. For brevity, 
we do not include the portions of code related to figure 
generation using the matplotlib framework [5J. 

4-1- API changes to dynamics solvers 

Here we demonstrate using the new Odedata class that 
is returned by the mcsolve, mesolve, brmesolve, and 
fmmesolve evolution solvers in QuTiP version 2.1. To bet- 
ter illustrate the API changes, we have recoded the two- 
qubit gate example from Ref. [T] Sec. (4.1) that is written 
using the older QuTiP 1.x API. The sections of the script 
featuring the new API are indicated below: 

from qutip import * 



g = 1.0 
gl = 0.75 
g2 = 0.05 
n_th =0.75 
T = pi/(4*g) 



2 * 



pi # coupling strength 

# relaxation rate 

# dephasing rate 

# bath avg. thermal excitations 

# gate period 



# construct Hamiltonian 

H = g * (tensor (sigmaxO , sigmaxQ) + 
tensor (sigmay () , sigmayO)) 

# construct inital state 

psiO = tensor (basis (2 , 1) , basis(2,0)) 

# construct collapse operators 
c_ops = [] 

## qubit 1 collapse operators 

sml = tensor (sigmamO , qeye(2)) 

szl = tensor (sigmazO , qeye(2)) 

c_ops . append (sqrt (gl * (l+n_th)) * sml) 

c_ops . append (sqrt (gl * n_th) * sml. dag O) 

c_ops . append(sqrt (g2) * szl) 

## qubit 2 collapse operators 

sm2 = tensor (qeye (2) , sigmamO) 

sz2 = tensor (qeye (2) , sigmazO) 

c_ops . append (sqrt (gl * (l+n_th)) * sm2) 

c_ops . append (sqrt (gl * n_th) * sm2.dag()) 

c_ops . append(sqrt (g2) * sz2) 

# evolve the dissipative system 
tlist = linspace(0, T, 100) 

medata = mesolve (H, psiO, tlist, c_ops, [] ) 
## NEW API CALL ## 

# extract density matrices from Odedata object 
rho_list = medata. states 

# get final density matrix for fidelity comparison 
rho_final = rho_list[-l] 

# calculate expectation values 

nl = expect (sml . dag() * sml, rho_list) 
n2 = expect (sm2 . dag () * sm2, rho_list) 

# calculate the ideal evolution 
medata_ideal = mesolve(H, psiO, tlist, [], [] ) 
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## NEW API CALL ## 

# extract states from Odedata object 
psi_list = medata_ideal . states 



0.6 r 



0.4 



# calculate expectation values 

nl_ideal = expect (sml . dagO * sml, psi_list) 
n2_ideal = expect (sm2 . dag () * sm2, psi_list) 

# get last ket vector for comparison 
psi_ideal = psi_list[-l] 

# output is ket since no collapse operators. 
rho_ideal = ket2dm(psi_ideal) 



< 0.2 



0.0 



O -0.2 



# calculate the fidelity of final states 
F = f idelity(rho_ideal, rho_final) 



-0.4 



4-2. Floquet modes of a driven two-level system 

Following the example in Ref. [7] , here we calculate the 
quasienergies for the time-dependent Floquet basis vectors 
of a sinusoidally driven two-level system [8] with Hamilto- 
nian 

A E 

H = —a z + — cos (ut) a x , (1) 

where A is the qubit energy splitting and to is the driving 
frequency, for different values of the driving amplitude E. 
The results of the simulation are presented in Fig. [I] . 

from qutip import * 

delta =1.0*2* pi # bare qubit sigma_z coefficient 

omega =8. 0*2* pi # driving frequency 

T = (2*pi) /omega # driving period 

# vector of driving amplitudes 

E_vec = linspace(0.0, 12.0, 100) * omega 

# generate spin operators 
sx = sigmaxO 

sz = sigmazO 

# create array for storing energy values 
q_energies = zeros ( (len(E_vec) , 2)) 

# define time-independent Hamiltonian term 
HO = delta/2.0 * sz 

args = { ' w ' : omega} 

# loop over driving amplitudes 
for idx, E in enumerate (E_vec) : 

# amplitude-dependent Hamiltonian term 
HI = E/2.0 * sx 

# H = HO + HI * cos(w * t) in 'list-string' format 
H = [HO, [HI, 'cos(w * t) ']] 

# find the Floquet modes 

f_modes, f_energies = f loquet_modes(H, T, args) 

# record quasi-energies 
q_energies [idx, :] = f_energies 



4-3. Floquet Evolution 

A driven system that is interacting with its environment 
is not necessarily well described by the standard Lindblad 
master equation as its dissipation process could be time- 
dependent due to the driving. In such cases a rigorous 
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Figure 1: (Color) Quasienergies corresponding to the two Floquet 
basis states of a driven two-level system as the driving strength is 
increased. Here the quasienergy and driving amplitude are expressed 
in units of the qubit energy splitting and driving frequency (h = 1), 
respectively. In this simulation, A = 1 X 2-7r and lu = 8 X 2ir. 



approach would be to take the driving into account when 
deriving the master equation. This can be done in many 
different ways, but one common approach is to derive the 
master equation in the Floquet basis, the Floquet-Markov 
master equation In QuTiP, this Floquet-Markov mas- 
ter equation is implemented in the fmmesolve function. 
As this approach is for time-dependent systems, here we 
model a sinusoidally driven qubit with Hamiltonian 



H 



A 



- a z — Aa x sin ujt 



(2) 



where A and e are the coupling and energy splitting con- 
stants, while A and lu are the driving strength and fre- 
quency, respectively. In addition, we define the spectral 
density of the environmental noise to be Ohmic. In Fig. [2] 
we plot the occupation probability of the qubit for both 
the Lindblad and Floquet-Markov master equations as a 
function of time. 

from qutip import * 



gammal = 0.05 
= 0.0 
= 0.0 * 2 
= 1.0 * 2 
= 0.1 * 2 
= 1.0 * 2 



delta 

epsO 

A 



T 

psiO 



# relaxation rate 

# dephasing rate 
pi # qubit sigma_x coefficient 
pi # qubit sigma_z coefficient 
pi # driving amplitude 
pi # driving frequency 

2*pi / w # driving period 

basis(2,0) # initial state 



tlist = linspaceCO, 25.0, 250) 

def J_cb(omega) : 

1111,1 Noise spectral density """ 
return 0.5 * gammal * omega/ (2*pi) 

# Hamiltonian in list-string format 
args = { ' w ' : w} 
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HO = - delta/2.0 * sigmaxO - epsO/2.0 * sigmazQ 

HI = - A * sigmaxO 

H = [HO, [HI, 'sin(w * t) ']] 



# Lindblad equation with time-dependent Hamiltonian 
# 

c_ops = [sqrt (gammal) * sigmaxO, 
sqrt(gamma2) * sigmazO] 
p_ex_me = mesolve(H, psiO, tlist, c_ops, 

[num(2)] , args=args) . expect [0] 



# Floquet-Markov master equation dynamics 
# 

rhs_clear() # clears previous time -dependent Hamiltonian 

# find initial Floquet modes and quasienergies 
f_modes_0, f_energies = f loquet_modes(H, T, args , False) 

# precalculate Floquet modes for the first driving period 
f _modes_table = f loquet_modes_table (f _modes_0 , f_energies, 

linspace(0, T, 500+1), H, T, args) 

# solve the Floquet-Markov master equation 
rho_list = fmmesolveCH, psiO, tlist, [sigmaxO] 

[] , [J_cb] , T, args). states 

# calculate expectation values in the computational basis 
p_ex_fmme = zeros (shape (p_ex_me) ) 

for idx, t in enumerate (tlist) : 

f_modes_t = f loquet_modes_t_lookup(f _modes_table, t, T) 
p_ex_f mme [idx] = expect (num(2) , 

rho_list [idx] . transf orm(f _modes_t , False)) 



— Lindblad 

— Floquet-Markov 




trace and positivity preserving), but it does not provide a 
connection to any underlying microscopic physical model. 
However, a microscopic model can in some cases be ad- 
vantagous, as for example in systems with varying energy 
biases and eigenstates that couple to an environment in 
some well-defined manner, through a physically motivated 
system-environment interaction operator that can be re- 
lated to a noise-power spectrum. The Bloch-Redfield for- 
malism is one such approach to derive a master equation 
from the underlying microscopic physics of the system- 
bath coupling. To highlight the differences inherent in 
these two approaches, in Fig. [3] we plot the expectation 
values for the spin-operators of the first qubit in a coupled 
qubit system given by the Hamiltonian 



H = LU 1 



cosflio-W + sin 0i ct^ 



(3) 



+ lo 2 cos0 2 al 2) +sin0 2 al 2) +ga^ai^ 



- w 2 

with the initial state 
5 



1^2) 



= |V>i) |^2) with 
[0.8 |0) (1) +(1-0.8)|1) (1) ; 
[(1 - 0.8) |0> (2) +0.8|1) (2) ; 



where the subscripts indicate which qubit the state belongs 
to. In Eq. pj), g is the qubit coupling, lo\ and u; 2 are 
the qubit frequencies, and finally 0i and 2 represent the 
angles of each qubit with respect to the o~ z direction. In 
this example, the qubit environments in the Bloch-Redfield 
simulation are assumed to have Ohmic spectrum. The 
code for the corresponding simulation is given below: 



10 15 
Time 



from qutip import * 

w = array ([1.0, 1.0]) *2*pi 
theta = array ( [0.025,0.0] )*2*pi 
gammal = [0.25, 0.35] 
gamma2 = [0.0, 0.0] 
g =0.1*2* pi 

# initial state 
a = 0.8 

psil = (a*basis(2,0) + (l-a)*basis(2,l)) .unitO 
psi2 = ((l-a)*basis(2,0)+a*basis(2,l)) .unitO 
psiO = tensor (psil, psi2) 

# times at which to evaluate expectation values 
tlist = linspace(0, 15, 500) 



# qubit angular frequency 

# angle from sigma_z axis 

# qubit relaxation rate 

# qubit dephasing rate 

# coupling strength 



Figure 2: (Color) Occupation probability of a sinusoidally driven 
qubit, initially in its ground state, under both Lindblad and Floquet- 
Markov master equation evolution, where the qubit parameters are 
A = 0, e = 1.0 X 27r, and the relaxation and dephasing rates are given 
by 71 = 0.05, 72 = 0. The driving term has amplitude A = 0.1 X 2tt 
and frequency ui = 1.0 X 2ir. Here, the spectral noise density of the 
environment is assumed to be Ohmic. 



4-4- Bloch-Redfield master equation 

The Lindblad master equation is constructed so that it 
describes a physical evolution of the density matrix (i.e., 



# operators for qubit 
sxl = tensor (sigmaxO 
syl = tensor (sigmayO 
szl = tensor (sigmazO 
sml = tensor (sigmamO 

# operators for qubit 
sx2 = tensor(qeye(2) , 
sy2 = tensor(qeye(2) , 
sz2 = tensor (qeye (2) , 
sm2 = tensor (qeye (2) , 

# Hamiltonian 

# qubit 1 

H = w[0] * (cos (theta [0] ) 

# qubit 2 

H += w[l] * (cos (theta [1] ) 



qeye (2)) 
qeye (2)) 
qeye (2)) 
qeye (2)) 

2 

sigmaxO ) 
sigmayO) 
sigmazO) 
sigmamO ) 



* szl + 



* sz2 + 



sin(theta[0] ) * sxl) 
sin(theta[l] ) * sx2) 
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# interaction term 
H += g * sxl * sx2 

# 

# Lindblad master equation 
# 

c_ops = [] 

c_ops . append(sqrt (gammal [0] ) * sml) 
c_ops . append(sqrt (gammal [1] ) * sm2) 

lme_results = mesolve(H, psiO, tlist, c_ops, 
[sxl, syl, szl]) 

# 

# Bloch-Redf ield master equation 
# 

def ohmic_spectruml (w) : 
if w == 0.0: 

# dephasing inducing noise 
return 0.5 * gamma2 [0] 

else : 

# relaxation inducing noise 

return 0.5 * gammal [0] *w/(2*pi) * (w > 0.0) 

def ohmic_spectrum2(w) : 
if w == 0.0: 

# dephasing inducing noise 
return 0.5 * gamma2[l] 

else : 

# relaxation inducing noise 

return 0.5 * gammal [1] *w/(2*pi) * (w > 0.0) 



1.0 




5 10 15 20 



brme_results = brmesolve(H, psiO, tlist, [sxl, sx2] , 
[sxl, syl, szl], [ohmic_spectruml , 
ohmic_spectrum2] ) 



4-5. Quantum process tomography 

To demonstrate quantum process tomography, here we 
simulate the effects of relaxation and dephasing on the 
two-qubit iSWAP gate [TO] when the qubits are coupled 
to a thermal environment with on average (n) = 1.5 ex- 
citations. The x-matrix obtained from QPT contains all 
the information about the dynamics of this open quan- 
tum system. In Fig. [4] we plot the x-niatrix for both the 
dissipative and corresponding ideal (unitary) iSWAP gate 
dynamics. 

from qutip import * 

g=1.0*2*pi # coupling strength 

gl = 0.75 # relaxation rate 

g2 = 0.25 # dephasing rate 

n_th =1.5 # bath avg. thermal excitations 

T = pi/(4*g) # gate period 

H = g* (tensor (sigmaxO , sigmaxO ) +tensor (sigmay () , sigmayO ) ) 
psiO = tensor (basis(2 , 1) , basis(2,0)) 

c_ops = [] 

# qubit 1 collapse operators 
sml = tensor (sigmamO , qeye(2)) 
szl = tensor (sigmaz () , qeye(2)) 

c_ops . append(sqrt (gl * (l+n_th)) * sml) 
c_ops . append(sqrt (gl * n_th) * sml.dagO) 
c_ops . append(sqrt (g2) * szl) 

# qubit 2 collapse operators 
sm2 = tensor (qeye(2) , sigmamO) 



Figure 3: (Color) Expectation values for the spin operators of qubit 
#1, for both the Lindblad master equation (top) and the Bloch- 
Redfield master equation (bottom), where the qubit environment in 
the latter case is assumed to have an Ohmic spectrum. Here, lui = 
UJ2 = 1 .Ox 'It,, 8\ = 0.025 X 2-7T, 82 = 0, and the qubit relaxation terms 
for qubit #1 and qubit #2 are give by 7^ = 0.25 and 7j 2 ' = 0.35, 
respectively. The qubit coupling is g = 0.05 X 2-7T. In this simulation, 
the dephasing terms are assumed to be zero. 



sz2 = tensor(qeye(2) , sigmazO) 
c_ops . append (sqrt (gl * (l+n_th)) * sm2) 
c_ops . append (sqrt (gl * n_th) * sm2.dag()) 
c_ops . append(sqrt (g2) * sz2) 

# process tomography basis 

op_basis = [[qeye(2), sigmaxO, sigmayO , sigmazO]] * 2 
op_label = [["i", "x", "y", "z"]] * 2 

# dissipative gate 

U_diss = propagator (H, T, c_ops) 
chi = qpt(U_diss, op_basis) 
qpt_plot_combined(chi , op_label) 

# ideal gate 

U_psi = (-lj * H * T).expm() 
U_ideal = spre(U_psi) * spost (U_psi . dag() ) 
chi = qpt(U_ideal, op_basis) 
qpt_plot_combined(chi , op_label) 



4-6. Exporting QuTiP data 

Finally, we demonstrate the exporting of data gener- 
ated in QuTiP to an external plotting program using the 
f ile_data_store and f ile_data_read to save and load 



G 




Figure 4: (Color) (a) The QPT process x-matrix for the dissipative 
iSWAP gate between two-qubits. Here the color indicates the phase 
of each matrix element. The qubit-qubit coupling strength is g = 
1.0 X 2-7T, whereas the relaxation and dephasing rates are gl = 0.75 
and g2 = 0.25, respectively, (b) The ideal iSWAP gate x _ma trix 
when the qubit dissipation and dephasing are not present. 



the data, respectively. To keep the example completely 
in Python, we have chosen to use Mayavi [TT] to plot a 
Wigner function in Fig. [5] generated in QuTiP correspond- 
ing to the state |^) = \a) + |— a) +\<fi) that is a Schrodinger 
cat state consisting of two coherent states with complex 
amplitude a, with an additional random ket vector |</>) 
created using QuTiP's random state generator: 

from qutip import * 

# Number of basis states 
N = 20 

# amplitude of coherent states 
alpha = 2.0 + 2.0j 

# define ladder operators 
a = destroy(N) 

# define displacement operators 
Dl = displace (N, alpha) 

D2 = displace (N, -alpha) 

# create superposition of coherent states + random ket 
psi = (Dl + D2) * basis(N.O) + 0.5 * rand_ket(N) 

psi = psi.unitO # normalize 

# calculate Wigner function 
xvec = linspace(-5, 5, 500) 
yvec = xvec 



W = wigner (psi, xvec, yvec) 

## new function calls ## 

# store Wigner function to file 

f ile_data_store ( "wigner .dat" , W, numtype= ' real ' ) 

# load input data from file 
input_data = file_data_readCwigner.dat') 

# plot using mayavi 

from mayavi. mlab import * 
X,Y = meshgrid(xvec , yvec) 

surf (xvec, yvec, input_data, warp_scale=' auto ' ) 

view(0, 45) 

showQ 




Figure 5: (Color) Wigner function for the state = \a) + |— a) + 
\<j>) where a = 2 + 2j is the coherent state amplitude, and \<f>) is a 
randomly generated state vector. This plot is generated in Mayavi 
using data from the QuTiP wigner function. 
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Appendix A. New functions in QuTiP 2.1 
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Quantum Object Methods 
conj 

eigenenergies 
groundstate 

matrix.element 

tidyup 

trans 

Bloch-Redfield Functions 
bloch jredf ield_solve 
bloch_redf ield_tensor 

Floquet Functions 

floquet _modes 

floquet _modes_t 

f loquet_modes_table 

floquet _modes_t _lookup 

f loquet_states 

floquet _states_t 

f loquet_state_decomposition 

floquet _wavef unction 

f loquet_wavef unction.t 

Evolution Solvers 

brmesolve 

f mmesolve 

Correlation Functions 

correlation 

correlation.ss 

Entropy /Entanglement Functions 
concurrence 
entropy .conditional 
entropyjnutual 

Quantum Process Tomography 

qpt 

qpt_plot 

qpt_plot -Combined 

Random State/ Operator Generation 

rand_dm 

randjierm 

rand_ket 

rand_unitary 



Conjugate of quantum object. 

Calculates the eigenvalues (eigenenergies if the operator is a Hamiltonian) for a quantum operator. 
Returns the eigenvalue and eigenstate corresponding to the ground state of the quantum operator. 

Gives the matrix element Q n m = (ipn |0| fpm) for the given operator and states tp„ and ij) m . 
Removes the small elements from a quantum object. 
Transpose of a quantum object. 



Evolve the ODEs defined by the Bloch-Redfield tensor. 

Bloch-Redfield tensor for a set of system-bath operators and corresponding spectral functions. 



Calculate the initial Floquet modes given a periodic time-dependent Hamiltonian. 

Calculate the Floquet modes at a time t. 

Calculate a table of Floquet modes for an interval of times. 

Look up the Floquet modes at an arbitrary time t given a pre-computcd Floquet-mode table. 

Calculate the initial Floquet states given a set of Floquet modes. 

Calculate the Floquet states for an arbitrary time t given a set of Floquet modes. 

Decompose an arbitrary state vector in the basis of the given Floquet modes. 

Calculate the initial wavefunction given a Floquet-state decomposition and Floquet modes. 

Calculate the wavefunction for an arbitrary time t given a Floquet-state decomposition and modes. 



Bloch-Redfield master equation solver. 
Floquet-Markov master equation solver. 



Transient two-time correlation function. 
Steady-state two-time correlation function. 



Calculate the concurrence entanglement measure for a two-qubit state. 

The conditional entropy S(A\B) = S(A,B) — S(B) of a selected density matrix component. 
Mutual information S(A : B) between selection components of a system density matrix. 



Quantum process tomography x-matrix for a given (possibly non-unitary) transformation matrix U. 
Visualize the quantum process tomography x-matrix. Plot the real and imaginary parts separately. 
X-matrix plot with height and color corresponding to the absolute value and phase, respectively. 



Random N X N density matrix. 
Random N X N Hermitian operator. 
Random iV X 1 state (ket) vector. 
Random N X N Unitary operator. 



Gates 
iswap 
sqrtiswap 

Utility Functions 
f ile_data_read 
f ile_data_store 
qload 
qsave 

rhs_generate 
rhs_clear 



Quantum object representing the iSWAP gate. 

Quantum object representing the square root iSWAP gate. 



Retrieves an array of data from the requested file 
Stores a matrix of data to a file to be read by an external program. 
Loads quantum object or array data contained from given file. 
Saves any quantum object or array data to the specified file. 

Pre-compiles the Cython code for time-dependent mesolve problems run inside a parf or loop. 
Clears string-format time-dependent Hamiltonian data. 



Table A.l: List of new user-accessible functions available in QuTiP 2.1. Additional information about each function may be obtained by 
calling help(function_name) from the Python command line, or at the QuTiP website [2]. 
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